Exact solution of the one-dimensional deterministic Fixed-Energy Sandpile 
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In reason of the strongly non-ergodic dynamical behavior, universality properties of deterministic 
Fixed-Energy Sandpiles are still an open and debated issue. We investigate the one-dimensional 
model, whose microscopical dynamics can be solved exactly, and provide a deeper understanding of 
the origin of the non-ergodicity. By means of exact arguments, we prove the occurrence of orbits 
of well-defined periods and their dependence on the conserved energy density. Further statistical 
estimates of the size of the attraction's basins of the different periodic orbits lead to a complete 
characterization of the activity vs. energy density phase diagram in the limit of large system's size. 



Sandpile models have been introduced by Bak, Tang 
and Wiesenfcld (BTW) as simple non-equilibrium sys- 
tems exhibiting self-organized criticality (SOC), a con- 
cept that provided the theoretical framework for the 
study of a large number of natural phenomena whose dy- 
namics self-organizes into a critical state in which long- 
range correlations and scaling laws can be detected Q, |i[ . 
More recently, the occurrence of SOC has been related 
to the presence of underlying absorbing state phase tran- 
sitions (APT) governing the critical behavior of the sys- 
tem 3]. The relationship between sandpiles and non- 
equilibrium APTs is emphasized by a class of models with 
locally conserved dynamics, the Fixed-Energy Sandpiles 
(FES) |^. Finite size scaling techniques led to the com- 
plete characterization of the universality class for driven- 
dissipative sandpile models with stochastic dynamics || , 
but failed in determining the critical properties of the 
deterministic BTW model, that have been partially re- 
covered by means of multiscaling analysis |6(. Stochas- 
tic FESs seem to belong to a different universality class 
with respect to directed percolation 0, Q, S @i while 
the deterministic BTW dynamics shows very strong non- 
ergodic effects, that cannot be understood using a purely 
statistical mechanics approach 0, • In dissipative sys- 
tems, the non-ergodicity of the BTW rule is related to the 
existence of recurrent configurations and to the abelian 
group structure of the configuration space uncovered by 
Dhar |lll | , and emerges when the process of grains addi- 
tion and dissipation is not completely random. In such 
cases, the configuration space breaks into separed regions 
in which the dynamics is gov erned by cyclic groups and 
turns out to be periodic 12]. In the limit of zero dis- 
sipation and addition, the group structure breaks down, 
but the system becomes similar to a FES, in which the 
non-ergodic behavior is recovered as the result of a self- 
sustained activity that eventually enters a periodic orbit. 

In this Letter, we report the exact solution of the one- 
dimensional deterministic Fixed-Energy Sandpile (ld- 
DFES), providing an exhaustive description of the micro- 
scopic dynamics and its steady states, and determining 
the system's macroscopic behavior in the limit of large 
size. Moreover, our explanation of the microscopic ori- 



gin for non-ergodicity in BTW FESs elucidates the re- 
lation, already suggested in Ref.^(J> with mode-locking 
phenomena in non-linear automata, i.e. when the sys- 
tem's dynamics blocks into orbits of fixed periodicity 
even though some external parameter (e.g. the energy) 
is continuously changed in a given interval (see, for in- 
stance, Ref.[l3J and references therein). 

Let us consider a ring of N sites, each one endowed 
with an energy z%, assuming non negative integer val- 
ues. Whenever the energy of a site equals or exceeds the 
threshold value z th — 2, the site becomes active and un- 
dergoes a toppling event, i.e. it looses 2 units of energy, 
that are equally redistributed among its 2 neighbors. On 
the contrary, if z\ < 2, the site i is stable. The system 
evolves with parallel dynamics, therefore the update rule 
can be written as follows: 

Zi(t + 1) - Zi(t) - 26[zi{t) - 2] + J2 %(*) - 2 ] (!) 

for i = 1,2, . . . ,N, where 6(x) = 1(0) if x > (x < 0) 
and Vt is the neighborhood of the site i. When the total 
energy E = J^. zi is sufficiently large, the conservation 
imposed by periodic boundary conditions prevents the 
system to relax in a configuration composed of all stable 
sites (absorbing state). Since the configuration space is 
finite and the dynamics is deterministic, after a transient, 
the system enters a periodic orbit (steady active state). 
We will conventionally assume sites to have energy be- 
tween and 3 (relaxing this condition does not modify 
the physical properties of the system), initial configu- 
rations being any possible sequence of N symbols ran- 
domly chosen in the alphabet {0,1,2,3}. In addition 
to the invariance under the action of the finite group 
of cyclic permutations and reflections on the ring, the 
system admits another internal symmetry: the configu- 
rations are dynamically equivalent under the transforma- 
tion Zi — > z'i = 3 — Zi. This means that the whole orbit 
traced starting from a configuration Z' — {z'i} is known 
if we know that traced by the configuration Z = {zi}, 
and the two limit cycles have the same period. 

In the mean-field description 0, , the critical behav- 
ior of the system around the absorbing-active transition 
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FIG. 1: The activity diagram p(£) (A) and the behavior of 
the period T as function of the energy density £ (B) for a 
ld-DFES of size N = 20. Figures display all results, without 
averaging, obtained starting from 100 randomly chosen initial 
configurations at each energy value. 



point is deduced studying the activity phase diagram, the 
plot of the activity P = jj ®( z i ~~ 2) (the density of 
active sites) vs. the energy density £ — E/N. The sym- 
metry of the problem allows to restrict the analysis of 
the ld-DFES to the energy density range £ £ [0, 3/2] (or 
in C S [3/2,3]). 

Numerical results show that large systems converge to 
periodic orbits of well-defined periods T in the whole 
range of the energy density. Fig. ^ reports the diagrams 
p(C) and T(£) for a system of size N = 20 and 100 ran- 
dom initial conditions. To emphasize the non-ergodicity, 
the results of all runs are displayed, without comput- 
ing any average. In particular, for £ < 1 the system 
converges to a completely frozen configuration (period 
T = 1) with no active sites, while in the range 1 < £ < 2 
only orbits of period T = 2 are observed. As imposed by 
symmetries, £ > 2 corresponds to fixed-point configura- 
tions with all sites active. At the critical points £ = 1,2, 
orbits of period N seem to be statistically dominant for 
large systems, but other periods (e.g. T = 1,2) are ob- 
served with lower frequency. 

The reasons for such a particular phase-diagram can 
be understood exploiting a result obtained in the math- 
ematical literature, in which BTW sandpiles have been 
studied on generic graphs with the name of parallel chip 
firing games (CFG) |lj]. When the underlying graph 
is undirected, with neither sinks nor sources, the CFG 
conserves the total energy. To our knowledge, the only 
result determining the properties of the periodic orbits 
of CFGs has been proved by Bitar and Goles in the 
case of trees. Their theorem states that on a tree the 
steady states of the BTW dynamics are fixed points or 
cycles of period two. If the graph contains some loops, 
the theorem does not hold. 

Notwithstanding, the method used in Ref.|l5| can be 
exploited to study the ld-DFES. Let us consider a system 
in a periodic steady state at a time to, and a temporal 
window P to = [to, t + T — 1] that corresponds to the first 
period from to . If the system entered the periodic orbit O 
at time ij n , the whole temporal support [i,- n , oc) of O will 
be indicated with supp(O). We define Sj(i) as a binary 



variable that assumes value 1 (0) if the site i is active 
(stable) at the time t. The sequence of all these binary 
values forms the activity vector Si — {si(t)\t £ supp(O)} 
of a site i. Moreover, a set [t, t +p — 1] C supp(O) is a 
maximal active set of length p > 1 if Sj(i + r) = 1 for 
r = 0,1, . . . ,p and s(t — 1) = s(t + p + 1) =0. Similarly, 
a set [t,t + q — 1] C supp(O) is a maximal stable set of 
length q > 1 if s(t-l) = s(t+q+l) = 1 and Si(t + r) = 
for r = 0, 1, ... , q. 

According to Lemma 1 of Ref . [lij . the following state- 
ment holds on a generic graph: if [s — k, s] C supp(O) is 
a maximal active (stable) set for a site i, then there exists 
a neighbor j of i such that [s — k — 1, s — 1] C supp(O) 
is a maximal active (stable) set for j. We call V (W) 
the maximum number of consecutive l's (0's) in all the 
activity vectors along the period, i.e. the length p (q) of 
the largest maximal active (stable) set over all sites. The 
analysis can be restricted to < V < T (0 < W < T), 
the cases V = (W = T) and V = T (W = 0) corre- 
sponding to fixed point configurations in which all sites 
are active (stable) [l^ . Moreover, because of the inter- 
nal symmetry with respect to ( — 3/2, VF's properties 
in the interval £ G [0,3/2] are equivalent to those of V 
in the interval £ £ [3/2,3]. As a consequence, we limit 
our study to W and corresponding maximal stable sets 
in the energy range C S [0, 3/2]. 

Two cases, W = 1 and W > 1, should be distinguished. 
If W = 1, the activity vectors of all sites in the system 
are 2-periodic. Suppose, indeed, that a site i topples at 
a time t together with m < 2 of its neighbors: at the 
following time step t + 1 the site i does not topple, but 
the remaining 2 — m neighbors topple. Thus, having lost 
2 — m energy units during the update at time t and gained 
the same quantity at the following time step, the value 
of site i has period 2. This argument holds for all sites in 
the periodic state (W = 1 for all sites), that consequently 
has global period T = 2. 

In order to study the case W > 1, suppose the largest 
maximal stable set [t,t + W — 1] is at a site Iq. For the 
above Lemma there is a neighbor %\ of io whose largest 
maximal stable set is [t— l,t+W — 2]. In turn, site i\ has 
a neighbor with largest maximal stable set [t — 2, t + W — 
3]. However, W > 1 imposes that i 2 ^ io, because the 
Lemma implies Si 2 (t — 1) = while from the definition 
of maximal stable set we need Sj (t — 1) = 1. Proceeding 
step by step along the ring in the same direction, we 
reach the site ijv = *o, whose maximal stable set is [t — 
N, t + W — N — 1] . The properties of the system in io (and 
in every other site) at time t — N result to be the same 
as that at time t, thus we can conclude that during this 
process, the system performs one or more periodic cycles 
and returns in the starting configuration after exactly N 
temporal steps. In other words, for W > 1 (V > 1) the 
period T divides N. 

The rest of this Letter provides a deeper insight into 
the structure of the periodic orbits using a simple method 
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FIG. 2: Graphical sketch of the method used to find out the 
structure of periodic configurations. The case W = 2 is con- 
sidered. Top figure shows the correspondence between the bi- 
nary representation of the spatio-temporal dynamics and the 
evolution of the real configuration. An example of the proce- 
dure used to recover the real values assumed by a given site 
during the period is reported in the bottom series of draws. 
Site values can be univocally determined in W + 1 steps. At 
each step, the value of the site at a certain time is computed 
from a number of constraints on its (spatio-temporal) neigh- 
bors. 



that is sketched in Fig. [21 for W = 2. Suppose that the 
interval [t, t + W — 1] of length < W < N is the largest 
maximal stable set for a site i, then Si(t + s) = for 
s = 0, ... ,W-1, but also Sj(i-l) = 1 said s^t+W) = 1. 
In addition, Si-i(t — 1 + r) = for r = 0, . . . , W — 1, with 
Si-i{t — 2) = 1 and Si_i(i + W — 1) = 1; and similarly 
s l+1 (t+l + r) = for r = 0, W-l, with s i+1 (t) = 1 
and s i+ i(t + W + 1) = 1. These binary values can be 
drawn on a spatio-temporal grid as in Fig. [21 In order 
to determine the real values assumed by i in the interval 
[t,t + W — 1], we need to discute separately the cases 
W = 1, W = 2 and W > 2. 

The case W = 1 is simpler because all maximal stable 
sets have length 1, thus Si(t) = implies that at least one 
of the two neighbors is active. Drawing a spatio-temporal 
grid as in Fig. [21 it follows that spatial configurations 
can be composed of all possible combinations of two-site 
blocks of the type 20, 21, 30 and 31 (and those with 
opposite order 02, 12, 03 and 13). Note that the set of 
all possible combinations of these blocks completely fills 
the energy interval ( € [1, 2] of 2-periodic orbits; while 
at the energy density £ = 1 (£ = 2), the configurations 
of the unique 2-periodic orbit are -Z20 = {• • • 202020 . . . } 
and Z 02 = {. . . 020202 . . . } (Z 31 = {. . . 313131 . . . } and 
Z 13 = {...131313...}). 

When W > 1 (see Fig. |2J), the site i does not topple 
at time t + W — 1 but it has to topple at time t + W, 
after having received a single energy unit (from the site 
i - 1) at time t + W - 1, then Zi(t + W - 1) = 1 and, 
consequently, Zi(t + W) = 2. In particular, if W = 2, the 



sites i and i — 1 do not topple at time t + W — 2 = t, 
but i + 1 topples, providing i of the unique unit of energy 
that we find at time t + W — 1 = t+ 1. Hence, Zi(t) = 0, 
from which follows that Zi(t — 1) = 2. 

In the case W > 2, the same arguments show that 
Zi(t + W — 2) = 1 (none of the neighbors provides any 
grain). If we go backward along the maximal comple- 
mentary set, a site i assumes value 1 up to a time t + 1, 
then Zi (t) — (one of its neighbors topples at that time) . 
In summary, during the time interval [t — 1, t + W] with 
W > 1, the site under study assumes a well-determined 
sequence of values 2011 ... 12 with exactly W — 1 values 
equal to 1. 

Fig. [21 shows the existing relation between the values 
assumed by a site during these temporal intervals and 
those admitted in the spatial arrangement of the corre- 
sponding periodic configurations. If the maximal stable 
set of length W for i starts at time t, the uniqueness of the 
previous construction allows to completely determine the 
spatial block {z i - 1 (t+W-l)z t (t+W~l) . . .z i+w (t+W- 
1)} that will be of the form {211 . . . 102} with W-l sites 
at 1. Applying a spatio-temporal translation to our con- 
struction, the spatial block {zi-w(t) ■ ■ ■ Zi(t) Zi + i(t)} — 
{211 . . . 102} is determined as well. 

The remaining structure of the configuration at time 
t can be studied with a similar technique, provided that 
we erase the spatial block {zi_w (t) . . . Zi(t)} and consider 
a reduced system composed of N — W — 1 sites. Since 
the role of the erased block was actually only that of 
transporting an activity 'soliton' from the site i — W to 
the site i + 1, this operation does not alter the dynamics 
that maintains a periodic behavior. After the reduction, 
the above methods are applied on the system, with in 
principle a largest maximal stable set of different length 
W < W. At the end of this process, the system (at 
a time t in the periodic regime) will be decomposed in 
blocks of decreasing length (from W + 1 to 2) and struc- 
ture of the type 211... 0. The only allowed blocks of 
length 2 are those of the form 20. The fronts direction 
of motion is determined by the initial conditions, but the 
dynamics is invariant under spatial reflections. Hence, 
the analysis can be applied to a system that evolves in 
the opposite spatial direction to that considered here, 
producing periodic configurations of the same structure 
but inverse spatial order. Moreover, the analysis of max- 
imal active sets for < V < N leads to identical results 
(in the range 3/2 < £ < 3) with spatial blocks of the 
type 122 .. . 31 (and 13 . . . 221) instead of 211 .. . 02 (and 
20... 112). 

We have proved by construction that the case W > 1 
(V > 1) corresponds exactly to systems with an energy 
density ( = 1 (£ = 2). For all other values of only 
periods of length 1 and 2 are allowed. In the region ( < 1 
(and C > 2) the limit cycles are fixed point configurations 
with all stable (active) sites, but in the interval 1 < C < 2 
fixed points are forbidden and W = 1. Then, in this 
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FIG. 3: In panel A), the number ^,t(N) of configurations 
belonging to the attraction's basin of orbits of period T — 1 
(crosses), 2 (squares), and N (circles) is plotted as function 
of N. The exponential rate of growth is extimated by curve 
fitting (dashed line for /ii(JV) and dot-dashed line for /ijv(A)). 
The inset shows that the ratio /ii(A)//ijv(A) decreases with 
N. Panel B) reports the total number fi(JV) of configurations 
at £ = 1 (or £ = 2) as function of the system size N (circles). 
The fit (dashed line) shows an exponential growth as fi(JV) tx 
a N with a ~ 3.5. 



region the period of the limit cycles is T = 2. Such exact 
result provides a theoretical foundation to the empirical 
data collected by numerical simulations (see Fig.QJ. 

At the critical points ( = 1 (( = 2), the knowledge 
of the dynamics of building blocks allows also to com- 
pute the exact number of periodic configuration of period 
T as a combinatorial enumeration problem. In particu- 
lar, the number H a u(N) of configurations belonging to 
orbits with T > 1 equals the number of ways a num- 
bered ring of N sites can be filled with an ordered set 
of blocks of length k comprised between 2 and N [lif. 
The generating function for the combinations with blocks 

of length 2 < k < N is Q(z) = l/(l - £ 2 <fc<jv z k ) ■ 
Hence, IL a u(N) is obtained as function of the coefficients 
of the derivatives in z — of Q(z), i.e. H a u(N) — 

T,2<k<N (N- k) \ £"n-% Q(z = 0) - A N , where A N = 
if A is odd and An — 2 if N is even. The last term 
is introduced to compensate the double counting (due 
to the factor 2k) of configurations composed of only 20 
blocks. The number n^(A) of configurations of period 
T < A is given by the number of different ways of fill- 
ing the ring with identical blocks of length and period 
equal to T, i.e. T1 T (N) = U T (T) for 2 < T < N. For- 
mally, Hat (AO = n a „(AT) - £ 2 < T<Ar n T (A) (with T di- 
visor of N). Starting from n 2 (2) = n aH (2) = 2 and 
113(3) = II a /;(3) = 6, a recursive relation allows to com- 
pute all the other terms IL T {N) (2 < T < N). The 
explicit computation shows that, in the large N limit, 
the mass of orbits of period T = N grows faster than 
that of orbits of periods T < N . 

However, in order to establish which period occurs 
with higher frequency at C = 1 (C = 2), the mass of 
the whole attraction's basin of an orbit is necessary, not 
only the number of periodic configurations. We have 
measured numerically the exact number of configurations 
in the basins of attraction of orbits of different peri- 



ods. Fig. UJ-A shows results up to JV = 16. The size 
/ii (N) of the basin of attraction of the fixed-point grows 
exponentially with a rate m(N) / fJ,i(N — 1) ~ £1 with 
£1 ~ 3.25. A similar growing rate is observed for or- 
bits of period T = 2 (A even), while the other orbits of 
period T < A present smaller attraction's basins. The 
largest growing rate is that of orbits of period A, for 
which fi N (N)/fi N (N - 1) ~ &v with £ N ~ 3.5. This 
means that for N — 100 the probability of observing 
T = 1 or T = 2 is about 10~ 4 smaller than that of ob- 
serving T — A. The conjecture that A-periodic orbits 
have probability 1 in the large A^ limit is corroborated 
by the observation of very similar scaling laws for /za? ( A~ ) 
and the total number d(N) of possible configurations at 
energy E = N, whose behavior is displayed in Fig. 0-B 
(data are computed analytically using simple combina- 
torics similar to that used in Ref. [lj|)- Consequently, 
at the critical energies £ = 1,2, with probability 1 the 
system enters very long orbits characterized by a steady 
current of active solitons transported along the system. 
The activity can assume all (rational) values between 
and 1/2. 

In conclusion, traditional statistical mechanics cannot 
capture the complex behavior of deterministic FESs, even 
in one-dimensional systems. On the contrary, our solu- 
tion provides the full understanding of the mode-locking 
phenomena and the sharp activity transitions observed in 
the phase diagram. The arguments can be partially ex- 
tended to higher dimensions, but the analysis is messed 
up by the complexity of the solitonic motion. In the 
case of two dimensional lattices, the present analysis to- 
gether with symmetry arguments allow to obtain an ap- 
proximated phase diagram 18], whose structure repro- 
duces the devil's staircase observed in numerical simu- 
lations 0. We hope that these results could also rep- 
resent a kind of benchmark for the analytical and nu- 
merical study of other automata with conserved dynam- 
ics, in which spatio-temporal patterns are governed by 
the motion of solitonic fronts. The author is grateful to 
M. Casartelli and P. Vivo for useful comments. 
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